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Abstract 

We investigate the nonequilibrium roughening transition of a one-dimensional 
restricted solid-on-solid model by directly sampling the stationary probabil- 
ity density of a suitable order parameter as the surface adsorption rate varies. 
The shapes of the probability density histograms suggest a typical Ginzburg- 
Landau scenario for the phase transition of the model, and estimates of the 
"magnetic" exponent seem to confirm its mean-field critical behaviour. We 
also found that the flipping times between the metastable phases of the model 
scale exponentially with the system size, signaling the breaking of ergodicity 
in the thermodynamic limit. Incidentally, we discovered that a closely related 
model not considered before also displays a phase transition with the same 
critical behaviour as the original model. Our results support the usefulness of 
off-critical histogram techniques in the investigation of nonequilibrium phase 
transitions. We also briefly discuss in an appendix a good and simple pseudo- 
random number generator used in our simulations. 
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1 Introduction 



During the last three decades or so, it has been found that nonequiUbrium one-di- 
mensional growth models may exhibit roughening transitions [1, 2, 3, 4], in con- 
trast with the widespread lore — based on simple entropy arguments — according 
to which one-dimensional interfaces in thermal equilibrium are always rough [5]. 
This discovery helped to expose, by means of concrete examples, some fundamen- 
tal differences that exist between equilibrium and nonequilium systems, e. g., that 
the dynamic fluctuations of a nonequiUbrium stationary state are in general not 
equivalent to the thermal fluctuations of an equilibrium state. It seems that, even in 
the stationary state, in nonequiUbrium interacting statistical systems time is more 
like another dimension of the system than a mere source of dynamic fluctuations, 
although in general it does not scale like the spatial dimensions of the system.^^ 

In Alon et al. [3] , a one-dimensional growth model was introduced that, when 
the adsorption of adatoms by the surface is small, desorption of adatoms from the 
edges of rough droplets together with absence of desorption from smooth terraces 
provide local mechanisms that eliminate islands of the rough phase, thus stabilising 
the flat phase despite the ergodicity of the model. Local rules dissolving islands 
of a minority phase provide an efficient mechanism of stabilisation against noise 
and have had a major role in the development of the foundations of nonequiUbrium 
statistical physics [6, 7, 8, 9, 10]. Recently, other conditions that can possibly give 
rise to phase-transitions in one-dimensional systems with short-range interactions 
have also been investigated; we refer the reader to ref. [1 1] for an exposition. 

In this work we investigate the roughening transition of the restricted solid-on- 
solid (RSOS) growth model introduced in ref. [3] by means of direct measurements 
of the stationary probability density of a suitable order parameter. In this way we 
were able to locate the critical point of the model and to estimate its "magnetic" 
critical exponent straighforwardly by Monte Carlo simulations. This sort of ap- 
proach seems to have been introduced in ref. [12] and has been used since then in 
varied contexts [13, 14, 15]. The break of ergodicity of the model in the thermody- 
namic limit could also be established from the flipping times between its metastable 
phases, which were found to scale exponentially with the system size. 

^That time could possibly be taken as an additional dimension in a related equilibrium model 
seems to have been suggested first by R. L. Dobrushin in mid-1960s to I. I. Piatetski-Shapiro and 
coworkers, that then began numerical work on cellular automata that ultimately led to the "positive 
probabilities conjecture." We are indebted to Professor A. L. Toom (UFPE) for this remark. 
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This article goes as follows. In section 2 we present the model we are interested 
in and our approach through the stationary probability density. In section 3 we 
detail our Monte Carlo simulations and estimate the critical point and the critical 
exponent of the model. In section 4 we look at an extended model in which 
desorption from smooth terraces are allowed and answer whether this extended 
model displays a phase transition and spontaneous symmetry breaking. Section 5 
discusses the flipping times between the metastable phases of the model and its 
ergodicity. Finally, in section 6 we summarise our conclusions and indicate some 
directions for further investigation. 



2 The RSOS growth model and the stationary probability 
density 

Let /i^ € N be the height of a surface or interface at site ^ of a finite lattice A C Z 
with |A| = L sites and periodic boundary conditions. The surface evolves by at- 
tempting, sequentially and at randomly chosen sites, adsorption of an adatom — )• 
/i£ + 1 with probability q dt, and desorption of an adatom hi — )■ min{/i^_i, hf\ or 
/i£ —7- min{/i£, /if+i} each with probability ^(1 — dt. We now impose the RSOS 
condition — h(\ ^ 1 G A, which suggests the use of the link variables 
Q = hij^i — hi G { — 1,0,1}. Denoting by F^^ the rate at which the elementary 
process (a, b) — )• (c, d) occurs, we can describe the above growth model in the 
links representation by the set of rates 

TZ = \{l-q). T% = q, (la) 

ro-° = g, T% = ]^[l-q), (lb) 
r+- = 1 - g, rf_ = g, (Ic) 

V = q- (Id) 

This model is translation invariant and conserves the total charge Q'^ — Q~ , with 
each sector of total charge corresponding to a closed class of the stochastic process. 
Notice that, differently from the models investigated in refs. [12, 15, 16, 17], this 
model does not conserve and individually. The adsorption and desorption 
moves corresponding to the above rates are represented in figure 1. 

According to (1), as q increases creation of H — as well as annihilation of — h 
pairs increase while the remaining processes induce spatial segregation of charges. 
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Figure 1: Elementary transitions in the RSOS growth model defined by the rates (1). The 
moving adatoms are depicted in grey. 

This corresponds to an increase in adsorption and in tiie growth of islands, leading 
to rougher configurations. As q lowers, increased annihilation of H — pairs together 
with more symmetric diffusion of O's flatten the surface. It is known that a rough- 
ening transition occurs at qc 
this transition is given by 

Ml 

This non-conserved order parameter anticipates the interpretation of the rough- 
ening transition as the spontaneous break of the Zoo symmetry related with the 
invariance of the growth process under an arbitrary integer shift hi ^ hi + n in 
the heights, for while in the rough phase all heights are exploited evenly, in the flat 
phase the system spontaneously selects one level around which the heights fluctu- 
ate. We then expect M = limi_j.oo Ml to be finite in the flat phase while vanishing 
in the rough phase due to canceling fluctuations. 

In this work we focus on the evaluation of the stationary probability density 
Pl{M) of the order parameter Ml defined above. Besides estimating the criti- 
cal point and the main critical exponent asociated with the transition directly from 
Pl{M), we also estimate the first passage time through a certain rough configura- 
tion from an initially flat configuration — the flipping time — and show that it grows 
exponentially with the system size, indicating that in the thermodynamic limit the 
process breaks ergodicity, signaling a phase transition. 



0.188 [3, 4, 18]. An order parameter that captures 
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A caveat is due about the approach employed in this article. Although it is 
tempting to define a free energy-like functional Fl{M) = —\nPi{M) once 
Pl{M) has been determined, we will avoid doing this. Firstly, because it is not 
necessary, since we do not intend to strugle with convexity and maximality issues 
here. Secondly, because it is not clear if the stationary measure can be expressed 
in the thermodynamic limit as a weighted exponential of some local function when 
detailed balance does not hold — Toom's model provides a striking counterexample 
[7, 8, 9]. Thus, except for one mention to a AM^ scenario in section 3, we avoid 
free energy-like arguments in what follows. 

3 Roughening transition 

Our Monte Carlo simulations ran as follows. The interface is initialised as a com- 
pletely flat interface with all heights hi equal, i. e., with all q = 0. In this con- 
figuration, the pseudo-spins mi should be initialised all with the same sign. It is a 
matter of choice to pick +1 or —1 for the initial mi, because we can anchor /ii(0) 
at any integer we like. 

For a given value of q, the initial configuration is relaxed through Monte 
Carlo steps (MCSs), with one MCS equal to L sequential attempts of movement. 
Ml is then sampled every other MCS and the data accumulated as a normalised 
histogram Pl{M). We drew the pseudo-random numbers used in our simulations 
from a combined linear congruential and 3-shift generator that is very fast, pos- 
sesses a period larger than 10^^, and passes several tests of randomness (cf. ap- 
pendix A) [19]. 

We obtained Pl (M) by using relatively small lattices and sampling a relatively 
large number of times. In order to obtain symmetric histograms, we change the 
sign of the pseudo-spins mg every after we sample M^. This allows us to sample 
both sides of the histogram without having to wait exponentially long times for 
the system to shift between them. This procedure is equivalent to simulating two 
systems in parallel. 

Scanning through q gives the set of curves shown in figure 2. This figure depicts 
a typical Ginzburg-Landau scenario of a second-order phase transition. Since we 
are dealing with a finite system, we could define Pl{M) ~ exp[—FL{M)] with 

Fl{M) = ^i^L{q)M^ + \xLiq)M^ + 0[M^), 
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Figure 2: Logarithm of the stationary probability densities P]^{M) for L = 512 and q 
between 0.180 (lowermost curve) and 0.188 (uppermost curve). Each curve was obtained 
from 10^ samplings of M^. We also depict a quadratic fit {y ~ —71.7 x^) to the minima 
of the curves (open circles), that although not being very tight, anticipates that the critical 
exponent 6 should be close to 1 /2. 

with fiLiq) = I^l\ic -Q) + OKQc - qf] and Xiiq) > 0. This form of Pl(M) 
would predict a mean-field critical behaviour for the order parameter, M oc {qc — 
qY with = \- This prediction is supported by previous results in the literature: in 
ref. [3], 6 was evaluated as 0.55(5) (the numbers between parentheses indicate the 
uncertainty in the last digits of the data), while for a model of yeastlike growth of 
fungi colonies with parallel dynamics it has been found that 9 ~ 0.50, although in 
this case the phase transition cannot be associated with the spontaneous breaking 
of a symmetry [20]. Moreover, in a certain line in the phase diagram of a closely 
related one-dimensional next-nearest-neighbour asymmetric exclusion process it 
has been found that 9 = 0.54 ± 0.04 [21]. In ref. [4], however, the shghtly less 
accurate but significantly higher value 9 = 0.66 ± 0.06 was published, pushing the 
estimate to a value that does not fit within a pure Ginzburg-Landau scenario. 

We estimated qc and 9 by fitting the modes of the probability histograms 

to 

Ml{q) = A{q-qc)'. (3) 

We tried to use the expected values {\Ml\) = Y.M Wl\Pl[M) instead of the 
modes, but they turned out to be too spread to be useful. The precision in is 
bounded from below by 5M = 2/L, the width of the bin of the probability his- 
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togram; in practice the uncertainties are higher because the histograms are some- 
what flat and noisy at the top. 

The modes for L = 512 are displayed in figure 3. For these data, if we fix 
= 5 we obtain the best fit with A = 2.752 and Qc = 0.1876, with an i?^ = 0.985 
and an adjusted R"^ = 1 - ^{1 - R'^) = 0.983, and if we fix qc = 0.1876 and 
let A and 6 vary we obtain A = 2.733 and 9 = 0.499, with the same and R"^ 
as before up to the fifth decimal place. These two fits cannot be discerned in the 
scale of figure 3 and are jointly indicated there by a solid line. Our data for larger 
lattices are shown in figure 4. Best fits of the L = 2048 data to (3) are obtained 
with A = 2.489 and qc = 0.1882 for fixed 6* = ^, with i?^ = 0.985 and adjusted 

^ 1 _ 1|(1 - i?2) ^ 0.983, and with A = 2.903 and 6 = 0.527 for q^ set at 
0.1882, now with i?^ = 0.999 and R"^ = 0.999. These values of A, 9, and qc are 
precise to the places given. 

4 The model with non-null desorption from terraces 

In the model analysed so far, local smooth terraces 00 can only evolve to dents 
H — by adsorption of an adatom. If we let T^_^_ > 0, this adsorption move gains 
the competition of a desorption move in which local smooth terraces 00 lose their 
middle adatom and become local pits — h. What happens when r^_^ > 0? 

The model with non-null desorption rate is relevant in the study of wetting (see. 



6 



0.2 




0.1 




1024 



1536 



□ 



X 



□ 



2048 



6=1/2 



e = 0.527 



0.02 



0.0001 



0.001 



0.01 



Figure 4: Logarithmic plot of A/£ for 512 ^ L ^ 2048. Each point was obtained from 
a histogram built from 10^ samplings of Af^. The solid and dashed lines indicate the best 
fits to the L — 2048 data for two different values of 9 (see text). 

e. g., ref. [22] and the references there in), where an additional condition is imposed 
on the system, namely, that there is a wall (or floor) over which the adatoms land 
or depart. There cannot be desorption from the wall, and this is a very strong 
constraint. In our formulation this constraint does not exist. 

There are many ways to take a finite T^^j^. One option — the most natural, in 
some sense, since the rate then becomes equal to the other desorption rate involving 
a smooth terrace — is to take T'^^j^ = 1—q, such that the rates T^^ and T^_^_^ become 
'dual' to the rates Fqq" and T^^_ and the model still has a single parameter q. After 
scanning through the whole interval ^ g ^ 1 for a system of L = 1024 sites, 
we did not find any sign of bimodality for Pl{M), and we conclude that with 
r'^'^_^ = 1 — q the growth model is always rough. We will not provide a graph 
of Pl{M) here because they are pretty featureless, relatively narrow distributions 
centered around M = 0. That in this case the model is always rough is somewhat 
expected, because the lower the adsorption of adatoms through the Foq^ process is, 
the more the complementary process T^^j^ digs pits in the interface, leading to an 
interface without any sizable smooth terrace. 

Another possibility to introduce a positive T^^j^ is to modify the original model 
(1) a little more hardly by setting Fqq^ = 1 — g and then adding a r'^'^_^ = q. In 
messing this way with the original model, we must first ask whether the modified 
model with T^^ = 1 — p displays a phase transition at all. The answer is yes. 
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with Qc ~ 0.131 and the same critical behaviour as the original model. Details 
about this modified model are not our concern here, though. We then proceed to 
modifying the original model as described above, add the desorption process with 
rate r^_^ = qto it, and seek a roughening transition. The result is that we did not 
find a roughening transition in this case either. 

Finally, we can assign an independent transition rate p £ (0, 1] to T^^. We 
then end up with a two-parameter model. We explored this model for small values 
of p. We found that even the smallest T^^, as small as p = 0.001 in a lattice of 
L = 1024 sites, suffices to destroy the flat order. 

As we have anticipated in the introduction (section 1), absence of desorption 
from smooth terraces provides a local mechanism that stabilises the flat phase. In 
the original model, the only way a flat terrace can be eroded is through desorption 
from its boundaries. Once we allow desorption from the middle of terraces, they 
loose their relative stability, particularly through the chain of reactions -| — — )• 
00 — )• — h, and does not survive anymore. 

5 Flipping times and ergodicity break 

The ergodicity of a nonequilibrium interacting particle system can be characterised 
by the exponential divergence of the flipping times between the different stationary 
states of the system as it gets larger. This divergence can be understood as the result 
of the superextensive growth of "energy barriers" between the metastable phases 
of the system as it gets larger, with the system getting trapped deeper inside one 
phase, ultimately leading to a break of ergodicity in the infinite system limit. 

Here we define the flipping time TL{q) as the first time (in MCSs) it takes the 
initial flat surface with Ml = 1 to become a rough surface configuration with 
Ml ^ 0. We take Ml = as the threshold because when the system reaches a 
surface configuration with Ml ^ coming from a configuration with Ml > 
it has, in an intuitive sense, "reached the other side of the well." Remember that, 
in principle, within each sector of total charge — Q~ of the stochastic process 
defined by rates (1) all configurations are reachable and the system is ergodic. 

For the purpose of exploring a possible break of ergodicity in model (1), it 
sufices to verify that TL{q) grows faster than algebraically with L. We thus look 
for a functional form 

rL(g) ~ exp(/L(g)). (4) 
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Figure 5: Flipping times as function of q and L. (a) The empirical {cf. (4)] grows 

as the system gets deeper inside the nonergodic phase q < qc- In this figure 0.155 ^ g ^ 
0.170 and L = 128 (lower curve), 160, 192, 224, and 256 (upper curve), (b) In this figure, 
q = 0.155 (upper curve), 0.157, 0.160, 0.165, and 0.170 (lower curve). In both cases we 
clearly see the exponential divergence of TL{q) as q | and as L ^ oo (for q < qc). 

A nonergodic dynamics implies that /L(g) diverges as L t cx) with finite q < qc, 
while for an ergodic dynamics /l(^J') should remain bounded in L. We also expect 
that T^iq) — )> oo as g I 0. This approach to determine the ergodicity of interacting 
particle systems has been applied successfully in several different models [23, 24]. 

We initialise the flat surface with all q = and Ml = L, release the system 
and count time until Ml ^ 0. We then obtained TL{q) for each pair of L and q as 
an average over 1000 such hitting times for some 128 ^ L ^ 256 and 0.155 ^ 
q ^ 0.170. These relatively small L and qc — q allow us to investigate TL{q) 
without having to wait too much to observe the flips. Our results appear in figure 
5. This figure clearly displays the exponential behaviour of TL{q) in the rough 
phase of the model, either as L t oo with fixed q < qc or as q \, for any L. 
This is strong indication that the model is nonergodic in the rough phase, meaning 
that in this phase the system does not sweep through all possible configurations 
compatible with the given sector of total charge — Q~ , getting clogged in 
some neighbourhood of the initial configuration. We also see from figure 5 (and 
log-log plots confirm) that seems to grow in a slightly superlinear fashion 

with {qc — q), i. e., as ~ {qc — qY with z >^\. Our data do not allow us to 

estimate z, though. 
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6 Summary and conclusions 



We showed that the roughening transition of the one-dimensional RSOS growth 
model given by rates (1) occurs at qc = 0.1880 it 0.0004 with an exponent 9 = 
0.51 ±0.02 associated with the order parameter Ml given by (2). The uncertainties 
in these figures are conservative, because we did not perform a proper finite-size 
scaling analysis of the model. Our value of 9 agrees with a couple of estimates 
published for this and related models [3, 20, 21], while disagreeing with ref. [4]. 
We believe that our data provide reasonable evidence for a mean-field exponent 
9 = ^. Incidentally, we discovered in section (4) that model (1) with a modified 
transition rate = 1 — q also displays a phase transition around qc ~ 0.131 
with the same critical behaviour as the original model. 

An analysis of the flipping times between the metastable phases of the model 
in the coarse-grained "M^-representation" enabled us to determine that the rough 
phase {q < qc) of the model corresponds to a nonergodic phase of the stochas- 
tic process, such that the phase transition of the model can be understood as an 
ergodic-nonergodic transition that does not find a counterpart in one-dimensional 
equilibrium systems at finite temperatures. 

Our results show that it is possible to investigate nonequilibrium phase tran- 
sitions by sampling the probability distribution of a suitable order parameter in 
the nonequilibrium stationary state of the process. While histogram techniques are 
standard in the simulation of equilibrium systems, for nonequilibrium systems they 
provide a viable alternative or complement to direct time-dependent simulations. 

Finally, we would like to mention that the order parameter Ml defined in (2) 
is actually one member of a family of order parameters given by 



These order parameters can be understood as discrete Fourier transforms of the 
probability distribution of the heights, and are closely related with the string or- 
der parameter introduced in ref. [25] to capture the disordered fiat phase of two- 
dimensional equilibrium crystal surfaces and the spin liquid phase (a type of dilute 
long-range antiferromagnetic order) of one-dimensional quantum spin chains of 
integer spin. These parameters behave in the RSOS model like[4] 




L 



(5) 



(n) 



'n 



(6) 
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each with a different 6'„. It would be interesting to estimate these exponents anew 
using our techniques to verify wether they assume their mean field values or not. 
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A The 'Simple KISS' pseudo-random number generator 

The 'Simple KISS' pseudo-random number generator (S-KISS for short) was in- 
troduced by G. Marsaglia in a newsgroup discussion as a simplification of his (and 
A. Zaman's) 'KISS' (Keep It Simple, Stupid) class of generators, which combine 
linear congruential, 3-shift, and multiply-with-carry generators [19]. The S-KISS 
generator keeps only the LCG and the 3-shift parts of its parent generator. Both 
generators pass some stringent statistical tests of randomness collectively known 
as the Diehard battery of tests [19]. 

In 32-bit arithmetic, the S-KISS generator is implemented by the following 
piece of code, intended to be included in the preamble of your own main program: 

static unsigned int A=32310901, C=1013904223; 
static unsigned int X=3593214833, Y=925021637; 
#define F32 (1 . 0/OxFFFFFFFF) 

#define rnd (X=A*X+C, Y"=Y<<13, Y"=Y>>17, Y'~=Y<<5, F32*(X+Y)) 

To draw a pseudo-random number just call rnd; e.g., p = rnd attributes 
a pseudo-random number (purportedly) distributed uniformly in [0, 1) to p. This 
generator has a period about as large as 2^^ ~ 1.8 x 10^^. 

Multiplier A is taken amongst the best for 32-bit LCG [26], while the initial 
values for C, X, and Y are arbitrary as long as C is odd and Y is non-null. Our initial 
X and Y have an equal number of bits and 1, not in any obvious pattern, just in 
case. You can change your favorite 32-bit LCG for the one given above. The shifts 
in the 3-shift part of the PRNG (the numbers 13, 17, and 5) as well as their order 
are not arbitrary; for other possible triples consult ref. [19]. 

*Thus i gnoring Einstein's advice to avoid making simple things simpler. 
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Neither the KISS nor the S-KISS generators were designed for cryptographic 
applications. This relative limitation notwithstanding, we believe that their sim- 
plicity, quality, high throughput, and large enough periods make them well suited 
for scientific applications like our Monte Carlo simulations. In fact, in this re- 
gard it would be interesting to examine these pseudo-random number generators 
using physics-based tests to assess their reliability in statistical physics simulations 
[27, 28]. 
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